﻿#include <pagmo/problem.hpp>
#include <pagmo/types.hpp>
#include <pagmo/algorithm.hpp>
#include <pagmo/algorithms/de.hpp>
#include <pagmo/population.hpp>
#include <iostream>
#include <vector>
#include <numeric>
#include <cmath>
#include <algorithm>
#include <random>
#include <iomanip>
#include <pagmo/algorithms/cstrs_self_adaptive.hpp>
#include <pagmo/algorithms/mbh.hpp>
#include <pagmo/algorithms/sade.hpp>
using namespace pagmo;

class CuttingStockProblem {
public:
    CuttingStockProblem() {};

    CuttingStockProblem(const std::vector<double>& demands,
        const std::vector<double>& lengths,
        double stock_length)
        : demands_(demands), lengths_(lengths), stock_length_(stock_length) {

        // 使用列生成算法思路生成切割模式
        generate_cutting_patterns_with_column_generation();
    }

    // 使用列生成算法思路生成切割模式
    void generate_cutting_patterns_with_column_generation() {
        patterns_.clear();

        // 1. 生成初始模式：每种零件单独切割（尽可能多放）
        std::cout << "Generating single-item patterns..." << std::endl;
        const int d_int = demands_.size();
        for (size_t i = 0; i < demands_.size(); ++i) {
            std::vector<int> pattern(d_int,0);
            pattern[i] = 1;
            patterns_.push_back(pattern);
        }

        // 2. 创建包含所有零件实例的大序列（每种零件重复25次）
        std::vector<size_t> all_item_instances;
        for (size_t i = 0; i < demands_.size(); ++i) {
            int count = static_cast<int>(demands_[i]);
            for (int j = 0; j < count; ++j) {
                all_item_instances.push_back(i);
            }
        }

        std::cout << "Total item instances: " << all_item_instances.size() << std::endl;

        // 设置随机数生成器
        std::random_device rd;
        std::mt19937 g(rd());

        // 3. 使用随机打乱+首次适应算法生成组合模式
        std::cout << "\nGenerating combination patterns with first-fit..." << std::endl;
        const int num_combinations = 100;
        auto calc_leng = [this](const std::vector<int>& pattern_tube) {
            double len = 0.0;
            for (size_t i = 0; i < pattern_tube.size(); ++i) {
                len += pattern_tube[i] * demands_[i];
            }
            return len;
        };
        for (int attempt = 0; attempt < num_combinations; ++attempt) {
            // 随机打乱所有零件实例
            std::shuffle(all_item_instances.begin(), all_item_instances.end(), g);

            // 使用首次适应算法生成切割模式
            std::vector<std::vector<int>> tubes;
            for (size_t idx : all_item_instances) {
                double item_length = lengths_[idx];
                bool found = false;
                for(auto& tube:tubes)
                {
                    const auto remaining_length = stock_length_ - calc_leng(tube);
                    // 检查是否还能放入该零件
                    if (item_length <= remaining_length) {
                        // 只能放入1个（因为每个实例代表1个零件）
                        tube[idx] += 1;
                        found = true;
                    }
                }
                if (!found) {
                    std::vector<int> current_pattern(demands_.size(), 0);
                    current_pattern[idx] += 1;
                    tubes.push_back(current_pattern);
                }
            }

            for (auto& tube : tubes) {
                add_pattern_if_unique(tube);
            }

            if (patterns_.size() % 20 == 0 && patterns_.size() > 0) {
                std::cout << "Generated " << patterns_.size() << " patterns so far..." << std::endl;
            }
            if (patterns_.size() > 100) {
                break;
            }
        }

        std::cout << "Total patterns generated: " << patterns_.size() << std::endl;
        print_pattern_statistics();
    }

    // 检查模式是否有效
    bool is_valid_pattern(const std::vector<int>& pattern) const {
        double total_length = 0.0;
        bool has_items = false;

        for (size_t i = 0; i < pattern.size(); ++i) {
            if (pattern[i] > 0) {
                has_items = true;
                total_length += pattern[i] * lengths_[i];
            }
        }

        return has_items && total_length <= stock_length_;
    }

    // 添加唯一模式
    void add_pattern_if_unique(const std::vector<int>& pattern) {
        for (const auto& existing_pattern : patterns_) {
            if (existing_pattern == pattern) {
                return;
            }
        }
        patterns_.push_back(pattern);
    }

    // 输出模式统计信息
    void print_pattern_statistics() const {
        std::cout << "\nPattern statistics:" << std::endl;
        std::cout << "===================" << std::endl;

        for (size_t i = 0; i < std::min(size_t(10), patterns_.size()); ++i) {
            std::cout << "Pattern " << i + 1 << ": ";
            double total_length = 0.0;
            int item_count = 0;

            for (size_t j = 0; j < patterns_[i].size(); ++j) {
                if (patterns_[i][j] > 0) {
                    std::cout << patterns_[i][j] << "x" << lengths_[j] << "mm ";
                    total_length += patterns_[i][j] * lengths_[j];
                    item_count += patterns_[i][j];
                }
            }

            double utilization = (total_length / stock_length_) * 100;
            std::cout << "-> " << item_count << " items, " << std::fixed << std::setprecision(1)
                << utilization << "% utilization" << std::endl;
        }
    }

    // 获取问题的维度
    vector_double::size_type get_nx() const {
        return patterns_.size();
    }

    vector_double::size_type get_nobj() const {
        return 1;
    }

    vector_double::size_type get_nec() const {
        return demands_.size();
    }

    vector_double::size_type get_nic() const {
        return 0; // 暂时不需要不等式约束
    }

    vector_double fitness(const vector_double& x) const {
        double total_stocks_used = 0.0;
        for (const auto& val : x) {
            total_stocks_used += val;
        }

        // 等式约束：每种零件的数量必须精确满足
        std::vector<double> eq_constraints(demands_.size(), 0.0);
        for (size_t j = 0; j < demands_.size(); ++j) {
            double total_produced = 0.0;
            for (size_t i = 0; i < patterns_.size(); ++i) {
                total_produced += x[i] * patterns_[i][j];
            }
            eq_constraints[j] = total_produced - demands_[j];
        }

        vector_double f(1 + get_nec());
        f[0] = total_stocks_used;

        for (size_t i = 0; i < eq_constraints.size(); ++i) {
            f[1 + i] = eq_constraints[i];
        }

        return f;
    }

    std::pair<vector_double, vector_double> get_bounds() const {
        vector_double lb(get_nx(), 0.0);
        vector_double ub(get_nx(), 50.0);
        return { lb, ub };
    }

    std::string get_name() const {
        return "Cutting Stock Problem with Instance-Based Patterns";
    }

    void print_solution(const vector_double& x) const {
        std::cout << "\nDetailed Cutting Plan:" << std::endl;
        std::cout << "======================" << std::endl;

        double total_stocks = 0.0;
        double total_used_length = 0.0;
        std::vector<double> produced(demands_.size(), 0.0);

        for (size_t i = 0; i < x.size(); ++i) {
            if (x[i] > 0.01) {
                double pattern_length = 0.0;
                std::cout << "Pattern " << i + 1 << " (used " << x[i] << " times): ";

                for (size_t j = 0; j < patterns_[i].size(); ++j) {
                    if (patterns_[i][j] > 0) {
                        std::cout << patterns_[i][j] << "x" << lengths_[j] << "mm ";
                        pattern_length += patterns_[i][j] * lengths_[j];
                        produced[j] += x[i] * patterns_[i][j];
                    }
                }

                double utilization = (pattern_length / stock_length_) * 100;
                std::cout << "-> " << std::fixed << std::setprecision(1) << utilization << "% utilization" << std::endl;

                total_stocks += x[i];
                total_used_length += x[i] * pattern_length;
            }
        }

        // 输出需求满足情况
        std::cout << "\nDemand Satisfaction:" << std::endl;
        std::cout << "====================" << std::endl;
        for (size_t i = 0; i < demands_.size(); ++i) {
            std::cout << "Item " << i + 1 << " (" << lengths_[i] << "mm): Required " << demands_[i]
                << ", Produced " << produced[i] << std::endl;
        }

        std::cout << "\nSummary:" << std::endl;
        std::cout << "========" << std::endl;
        std::cout << "Total stocks used: " << total_stocks << std::endl;
        std::cout << "Total material used: " << total_used_length << "mm" << std::endl;
        std::cout << "Overall utilization: " << (total_used_length / (total_stocks * stock_length_)) * 100 << "%" << std::endl;
    }

private:
    std::vector<double> demands_;
    std::vector<double> lengths_;
    double stock_length_;
    std::vector<std::vector<int>> patterns_;
};

int csp() {
    std::vector<double> demands = { 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
                                   25, 25, 25, 25, 25, 25, 25, 25, 25, 25 };
    std::vector<double> lengths = { 120, 270, 330, 450, 520, 590, 660, 740, 980, 1010,
                                  1250, 1380, 1560, 1780, 2540, 3080, 3520, 3900, 4090, 4500 };
    double stock_length = 6000;

    CuttingStockProblem problem(demands, lengths, stock_length);
    pagmo::problem prob{ problem };
    // 替换算法部分

    //pagmo::algorithm algo{ pagmo::de(1000, 0.8, 0.9, 2, 1e-6, 1e-6) };
    pagmo::algorithm algo{ pagmo::cstrs_self_adaptive(1000) };
    // 或者
    // pagmo::algorithm algo{pagmo::sade(1000)};
    pagmo::population pop{ prob, 200 };

    std::cout << "Starting optimization..." << std::endl;

    for (int i = 0; i < 100; ++i) {
        pop = algo.evolve(pop);

        if ((i + 1) % 20 == 0) {
            auto best_f = pop.champion_f();
            double eq_violation = 0.0;
            for (size_t j = 1; j < 1 + prob.get_nec(); ++j) {
                eq_violation += std::abs(best_f[j]);
            }
            std::cout << "Generation " << (i + 1) << ": Stocks = " << best_f[0]
                << ", EQ violation = " << eq_violation << std::endl;
        }
    }

    auto best_solution = pop.champion_x();
    problem.print_solution(best_solution);

    return 0;
}